from cmath import log
import math
from pstats import Stats
import numpy as np
import pandas as pd
import scipy.optimize as optimize
import scipy.interpolate as spi
import random
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
import scipy.stats as stats
from scipy import interpolate
import statsmodels.api as sm
import warnings

#---------------------------------------定义故障分布函数

'''
利用三次样条插值法计算故障率
'''
def func_cdf_cub(x, a,b,c):   
    t = (a,b,c) 
    return interpolate.splev(x,t)

'''
计算高斯分布时的故障率
'''
def func_cdf_gau(x, miu, sigma):    
    return 1-stats.norm.cdf((x-miu)/sigma)

'''
计算对数高斯分布时的故障率
'''
def func_cdf_loggau(x, miu, sigma):    
    return 1-stats.norm.cdf((np.log(x)-miu)/sigma)

'''
计算威布尔分布时的故障率
'''
def func_weib(x, m, gamma, t0):
    try:
        x = np.array(x)
        return np.exp(-1 * (x - gamma) ** m / t0)
    except OverflowError:
        return "该部件故障率不符合威布尔分布，请重新选择"

'''
计算指数分布时的故障率
'''
def func_expon(x, b):
    return np.exp(-b / x)

'''
拟合指数分布
'''
def fit_expon(middle, numbers):
    popt, pcov = optimize.curve_fit(func_expon, middle, numbers)
    return popt[0]

'''
拟合高斯分布
'''
def fit_gau(middle, numbers):
    popt, pcov = optimize.curve_fit(func_cdf_gau, middle, numbers)
    return popt[0], popt[1]

'''
拟合对数高斯分布
'''
def fit_loggau(middle, numbers):
    popt, pcov = optimize.curve_fit(func_cdf_loggau, middle, numbers)
    return popt[0], popt[1]

'''
拟合威布尔分布
'''
def fit_weib(middle, numbers):
    loc = stats.exponweib.fit(numbers, floc=0, f0=1)
    return loc[1], loc[2], loc[3]

'''
三次样条插值法拟合
'''
def fit_cub(middle, numbers):
    t = interpolate.splrep(middle, numbers, k = 3)
    a = t[0]
    b = t[1]
    c = t[2]
    return a,b,c

'''
主函数
'''
def run(input_array, type = 5, ct = 0):
    input_array = np.sort(input_array) #原始数据排序
    min = input_array[0]#整体下界 
    max = input_array[-1]#整体上界
    unit_number = 25 #分块个数
    np.set_printoptions(precision = 3, suppress = True)
    unit_length = (max - min) / unit_number #单位长度
    upper_border = [] #区块上界
    lower_border = [] #区块下界
    middle = [] #区块平均值，用作横坐标
    lower = min 
    upper = min + unit_length
    numbers = [] #每个区块内的样本个数，用作纵坐标

    #————————————————————————————————确定每个区块的左右边界
    for i in range(0,unit_number):
        lower_border.append(lower)
        upper_border.append(upper)
        middle.append((lower + upper)/2)
        lower = upper
        upper += unit_length
        numbers.append(0)

#——————————————————————————————————————计算每个区块内的样本个数
    index = 0    
    for i in range(0, len(input_array)):
        while (index < unit_number - 1) & (input_array[i] > upper_border[index]):
            index += 1
        numbers[index] += 1

    for i in range(1, unit_number +1):
        numbers[i-1] = numbers[i-1]/len(input_array)

    index = 0
    for i in range (0, len(middle)):
        if numbers[index] == 0:
            del middle[index]
            del numbers[index]
            index = index - 1
        index += 1
    middle = np.array(middle)
    numbers = np.array(numbers)

#——————————————————————————————————————————————————————————分布拟合
    if type == 1:   
        #print({'指数分布参数：':(fit_expon(middle, numbers))})
        result = (fit_expon(middle, numbers))
        print({"code": 1, "msg": "指数分布参数", "params": fit_expon(middle, numbers).tolist(), "chart": {"x": x, "y": [round(i, 3) for i in func_expon(x, fit_expon(middle, numbers))]}, "gzl": func_expon(ct, result)}, end = "")

    elif type == 2:   #高斯分布
        a,b = fit_gau(middle, numbers)
        #print({'高斯分布参数：':a, '参数b：':b})
        print({"code": 1, "msg": "高斯分布参数", "params": list(fit_gau(middle, numbers)), "chart": {"x": x, "y": [round(i, 3) for i in func_cdf_gau(x, a, b)]}, "gzl": func_cdf_gau(ct, a, b)}, end = "")

    elif type == 3:   #对数高斯分布
        a,b = fit_loggau(middle, numbers)
        #print({'对数高斯分布参数：':a, '参数b：':b})
        print({"code": 1, "msg": "对数高斯分布参数", "params": list(fit_loggau(middle, numbers)), "chart": {"x": x, "y": [round(i, 3) for i in func_cdf_loggau(x, a, b)]}, "gzl": func_cdf_loggau(ct, a, b)}, end = "")

    elif type == 4:    #威布尔分布
        warnings.filterwarnings("ignore")
        a,b,c = fit_weib(middle, numbers)
        #print({'威布尔分布参数a：':a, '参数b：':b, '参数c：':c})
        print({"code": 1, "msg": "威布尔分布参数", "params": list(fit_weib(middle, numbers)), "chart": {"x": x, "y": [round(i, 3) for i in func_weib(x, a, b, c)]}, "gzl": func_weib(ct, a, b, c)}, end = "")

    else: #三次样条差值法
        a, b, c = fit_cub(middle, numbers)
        #print({'三次样条插值法参数a：':a, '参数b：':b, '参数c：':c})
        gzl = func_cdf_cub(ct, a, b, c)
        print({"code": 1, "msg": "三次样条插值法", "a": [round(i, 2) for i in a], "b": [round(i, 2) for i in b], "c": c, "result": [round(i, 2) for i in func_cdf_cub(x, a, b, c)], "chart": {"x": x, "y": [round(i, 3) for i in func_cdf_cub(x, a, b, c)]}, "gzl": float(gzl)}, end = "")
        '''
        三次样条插值法的输出为元组（a,b,c），若计算时请一并将a,b,c按照顺序传回程序中的func_cdf_cub（x,a,b,c）函数中获取计算结果
        '''

#————————————————————————————
#指数分布求剩余工作时间
def find_expon(reliability, b):
    return -b * log(reliability)

#——————————————————————————————
#高斯分布求剩余工作时间
def find_gau(reliability, miu, sigma):
    reliability = 0.5
    worktime = 0.1
    low = 0
    high = worktime
    while (func_cdf_gau(high,miu, sigma) > reliability):
        high = high *2
    middle = (high + low) / 2
    count = 0
    while (abs(func_cdf_gau(middle, miu, sigma)-reliability) > 0.01 and count < 1000):
        if (func_cdf_gau(middle, miu, sigma) < reliability):
            high = middle
        else:
            low = middle
        middle = (high + low) / 2
        count += 1
        print(middle)
    middle = round(middle, 4)
    return middle

#————————————————————————————————
#对数高斯分布求剩余工作时间
def find_loggau(reliability, miu, sigma):
    reliability = 0.5
    worktime = 0.1
    low = 0
    high = worktime
    while (func_cdf_loggau(high,miu, sigma) > reliability):
        high = high *2
    middle = (high + low) / 2
    count = 0
    while (abs(func_cdf_loggau(middle, miu, sigma)-reliability) > 0.01 and count < 1000):
        if (func_cdf_loggau(middle, miu, sigma) < reliability):
            high = middle
        else:
            low = middle
        middle = (high + low) / 2
        count += 1
        print(middle)
    middle = round(middle, 4)
    return middle

#——————————————————————
#威布尔分布求剩余工作时间
def find_weib(reliability, m, gamma, t0):
    return gamma + pow((-t0*log(reliability), 1/m))

#——————————————————————————
#三次样条插值法求剩余工作时间
def find_cub(reliability, a,b,c):
    reliability = 0.5
    worktime = 0.1
    low = 0
    high = worktime
    while (func_cdf_cub(high,a,b,c) > reliability):
        high = high *2
    middle = (high + low) / 2
    count = 0
    while (abs(func_cdf_cub(middle, a,b,c)-reliability) > 0.01 and count < 1000):
        if (func_cdf_cub(middle, a,b,c) < reliability):
            high = middle
        else:
            low = middle
        middle = (high + low) / 2
        count += 1
        print(middle)
    middle = round(middle, 4)
    return middle

# 部件级态势感知算法交付
if __name__ == "__main__":
    #run([1,2,3,4], type = 3)
    #print(func_weib(5, 3446685051.4268165, 0, 0.24999999999552736))
    tp = 5
    x = [1,2,3,4]
    ct = 5.1
    #run([1,2,3,4], type = tp) # 第一个参数x是选择该部件的所有历史故障时间数据，只能选择一个部件，第二个参数type是计算类型，第三个参数ct是部件当前工作时间，直接部件故障率算出来
    run(x, type = tp, ct = ct)
    #print(func_expon(x,2))
    #print(func_weib(5, 3446685051.4268165, 0, 0.24999999999552736))
    
    '''
    type = 5 三次样条插值法
    '''
    #print(func_cdf_cub(5, [1.06, 1.06, 1.06, 1.06, 3.94, 3.94, 3.94, 3.94], [0.25, 0.25, 0.25, 0.25, 0.  , 0.  , 0.  , 0.  ], 3))

    '''
    type = 4 威布尔分布
    '''
    #print(func_weib(0.5, 3446685051.4268165, 0, 0.24999999999552736))

    '''
    type = 3 对数高斯分布
    '''
    #print(func_cdf_loggau(5, -2.436e+08, 3.611e+08))

    '''
    type = 2 高斯分布
    '''
    #print(func_cdf_gau(5, -2.436e+08, 3.611e+08))

    '''
    type = 1 指数分布
    '''
    #print(func_expon(5, 3.924))